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Abstract 

We consider the Klein-Gordon and sine-Gordon type equations with a point- 
hke potential, which describes the wave phenomenon in disordered media with 
a defect. The singular potential term yields a critical phenomenon-that is, 
the solution behavior around the critical parameter value bifurcates into two 
extreme cases. Finding such critical parameter values and the associated statis- 
tical quantities demands a large number of individual simulations with different 
parameter values. Pinpointing the critical value with arbitrary accuracy is even 
more challenging. In this work, we adopt the generalized polynomial chaos 
(gPC) method to determine the critical values and the mean solutions around 
such values. 

First, we consider the critical value associated with the strength of the singu- 
lar potential for the Klein-Gordon equation. We show the existence of a critical 
behavior with certain boundary conditions. Then we expand the solution in 
the random variable associated with the parameter. The obtained partial dif- 
ferential equations are solved using the Chebyshev collocation method. Due 
to the existence of the singularity, the Gibbs phenomenon appears in the solu- 
tion, yielding a slow convergence of the numerically computed critical value. To 
deal with the singularity, we adopt the consistent spectral collocation method. 
The gPC method, along with the consistent Chebyshev method, determines the 
critical value and the mean solution highly efficiently. 

We then consider the sine-Gordon equation, for which the critical value 
is associated with the initial velocity of the kink soliton solution. The critical 
behavior in this case is that the solution passes through (particle-pass) , is trapped 
by (particle-capture), or is reflected by (particle- reflection) the singular potential 
if the initial velocity of the soliton solution is greater than, equal to, or less than 
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the critical value, respectively. Due to the nonlinearity of the equation, we use 
the gPC mean value rather than reconstructing the solution to find the critical 
parameter. Numerical results show that the critical value can be determined 
efficiently and accurately by using the proposed method. The results are also 
compared with the results using the Monte-Carlo method. 
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1. Introduction 

Critical phenomena of partial differential equations (PDEs) associated with 
the criticalparameter can be found in various physical applications. For ex- 
ample, in |5| the gravitational collapse of the scalar field was considered. The 
critical phenomenon of the gravitational collapse was discovered; there exists a 
critical initial mass such that if the initial mass is less than the critical mass, the 
scalar field eventually disperses. Conversely, it forms a black hole if the value 
exceeds the critical mass. If the initial parameter is arbitrarily close to the crit- 
ical value, scale-invariant behavior of the solution occurs 5]. For PDEs from 
nonlinear optics, a similar phenomenon has been found and studied, e.g. the 
sine-Gordon equation or nonlinear Schrodinger equation with a point-like sin- 
gular potential term l^, H, [3 [lij . For the sine-Gordon equation, there exists 



a kink soliton solution; this yields a critical behavior once a self-interacting sin- 
gular potential is placed. This phenomenon is known as the wave phenomenon 
in disordered media with a defect. If the initial velocity of the kink soliton 
solution exceeds the critical velocity, the kink soliton eventually passes through 
the potential, which is known as the particle-pass. If the initial velocity is less 
than the critical velocity, the soliton is either captured by or reflected by the 
potential, known as the particle- capture or particle-reflection, respectively. Since 
the sine-Gordon equation is not integrable when the singular potential term is 
present, no exact solution is available; no exact value of the critical parameter is 
available either. There are many unknown properties of the solution behaviors 
in this case. For more detailed study, we refer the readers to 13, 12, 1^ 14. 15|. 



A similar phenomenon has been also found with nonlinear Schrodinger equa- 



tions 



26 . 30| . For example, in [sl [26| the singular potential term perturbs the 
soliton propagation. Similar phenomena of particle-pass, particle-capture, and 
particle-reflection were observed for some critical parameters. 

In this paper, we consider the critical phenomenon that is induced by the 
point-like singular source term in the PDEs, such as the sinc-Gordon equation 
mentioned above. In particular, we consider the Klein-Gordon and sine-Gordon 
equations to be perturbed by the singular potential term. For both cases, the 
solution behaviors are highly sensitive to the critical value of the parameter. 
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Determination of the critical value and the mean solution is an important task 
because the wave or soliton interactions with the local defect should be well 
understood for better optical communications. In either particle-capture or 
particle-reflection, the optical signal is not transmitted with the associated para- 
metric value if it is less than the critical value. The signal should be transmitted 
with some values of the parameter that exceed the critical value for the com- 
munication without blocking or trapping 27|. Finding the critical value is 
doable but challenging. 

There are two reasons. First, the singular potential is given by the Dirac 
(5-function, so the given PDEs are not well-defined. For the sine-Gordon equa- 
tion, no exact solution is known with the presence of the singular potential 
term, which means that no analytic form of the critical value is available-only 
the numerical or real laboratory experimental approach can render the criti- 
cal value measurable. The numerical determination of the critical parameter 
value relies on expensive Monte-Carlo(MC) simulations based on the bi-section 
method. Moreover, it is highly expensive to determine the critical value up to 
an arbitrary accuracy (for the sine-Gordon case, it is not even known whether 
the critical value is rational or irrational). Furthermore, the numerically ob- 
tained critical value depends on the numerical method, considered with given 
numerical settings such as the truncation number, domain size, final time etc. 



3l|. Due to the singularity in the potential term, the Gibbs phenomenon also 



appears in the finite truncated solution, resulting in a slow convergence. Thus 
when we talk about finding the critical value, we mean the critical value defined 
by the numerical method used, not the exact value from the PDEs. This will 
be explained in detail, using the Klein-Gordon equation in the paper. 

Second, although the singular potential term is used to mimic the effect of 
the local defect, such a model is not perfect. Associated parameters, such as the 
strength and the location of the singular potential terms, have a great degree 
of uncertainty. Thus, determining the critical value needs to be understood 
within the context of the uncertainty analysis. In this paper, as the first step 
toward thoroughly analyzing these two issues, we propose to use the generalized 
polynomial chaos (gPC) method ll|, |2J, [SSj for finding the critical value of 
both the Klein-Gordon and sine-Gordon equations and the associated statistical 
quantities such as mean solutions around the critical value. The gPC method 
has been successfully implemented to analyze the uncertainty quantification for 
problems such as fluid flows with uncertain initial and boundary conditions, 
sensitivity analysis, steady-state problems, etc. [37, 

In our proposed method, we consider the associated parameter that yields 
the critical phenomenon of the system as a random variable, and assume that 
we will try to find a critical parameter value with the given numerical setting. 
For the Klein-Gordon equation, the associated parameter is the strength or 
amplitude of the singular potential term; for the sine-Gordon equation, the as- 
sociated parameter is the initial velocity of the soliton solution in the presence 
of the singular potential term that has a fixed amplitude. The original equa- 
tions are then parameterized by the random variable and the solution of the 
parameterized equations is obtained by the gPC method. That is, the solution 
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is expanded by the orthogonal polynomials in the random space. The associ- 
ated orthogonal polynomials are determined by the probabilistic nature of the 
random variable, which is the critical parameter in our case. For example, if the 
random variable has the uniform probability density function (PDF), the associ- 
ated polynomial is the Legendre polynomial; if the PDF is normal distribution, 
the associated polynomial is the Hermite polynomial [iJ, [s^ ■ 

We first consider the Klein-Gordon equation with the singular potential and 
find some conditions for the critical phenomenon. The Klein-Gordon equation 
with the singular source term has been considered in various studies. For ex- 
ample, in [23|, the singular source term is considered in adopting the model, 
the meson theory, that the particle mass varies with space and time. In Q, two 
(5-functions are used as the singular potentials for Schrodinger equation and the 
corresponding eigenfunctions are derived using the Weyl-Titchmarsh-Kodaira 
(WTK) spectral theorem. In [13], the Klein-Gordon equation with the singu- 
lar defect was considered and the chaotic behavior of the solution is studied. 
However, not every singular potential term yields the critical phenomenon. We 
consider the self-interacting singular potential term and show the existence of 
the critical phenomenon with a certain type of boundary condition. That is, the 
energy norm increases, vanishes, or decreases if the parameter value exceeds, is 
equal to, or is less than the critical value. The critical value of the parameter is 
obtained analytically in this case, which makes it possible to do an error anal- 
ysis. For the critical value, we assume the uniform PDF and use the Legendre 
polynomial for the expansion. The resulting deterministic equations are solved 
using the Chebyshev collocation method with the consistent formulation for the 
Dirac (5-function [l^ • Once the expansion coefficients are obtained, the solution 
is reconstructed for various values of the parameter. Using the fact that the 
energy norm remains constant for the critical value when the time is large, the 
invariant point is found as the critical value. The mean solution is also found 
by the first gPC mode. We compare the gPC approach and the MC approach. 

Then we apply the same method for the sine-Gordon equation. The sine- 
Gordon equation with singular potential terms is also found in much of the 
existing literature. The singularly perturbed solution of the sine-Gordon equa- 
tion shows many interesting mathematical structures, such as the fractality, 



self-similarity, etc. [12|, [13|. Moreover, such an equation yields the critical phe- 
nomenon associated with the critical parameter, as explained above briefly. In 
our work, we use the initial velocity of the soliton solution as the random vari- 
able, and try to find the critical initial velocity around which the particle-pass, 
particle-capture, and particle-reflection occur. For the random variable, we con- 
sider both the Hermite and Legendre polynomials-that is, the normal and the 
uniform distributions as the PDF. Unlike the Klein-Gordon equation, the recon- 
struction of the solution using the gPC method is not obtained directly. This is 
because the receiving soliton signal has two extreme cases, i.e. transmitted or 
not-transmitted, and because of the nonlincarity. The extreme solution behav- 
ior renders the PDF as the two (5-function-likc PDFs. For the particle-capture, 
the amplitude of the tail at the domain exit is 27r; it vanishes eventually for 
the particle-pass. Such strong discontinuity in the solution behavior causes the 
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Gibbs phenomenon in the reconstruction. These two extreme cases are not 
separable in the reconstruction. In the paper, we show such non-separable re- 
construction for both the Hermite and Legendre cases. Thus instead of using 
the reconstruction by the gPC method, we use the gPC mean for the determi- 
nation of the critical value. In this case, the Legendre gPC method yields much 
faster convergence than the Hermite gPC method. Numerical results show that 
the Legendre gPC method determines the critical value very efficiently and ac- 
curately. The determined critical value agrees well with the value found by the 
many direct MC simulations. 

To the best of our knowledge, gPC analysis for the singularly perturbed 
PDEs by the Dirac (5-functions has not been thoroughly studied in existing 
literature. In particular, for the nonlinear optics equations with a singular 
defect in the disordered media, the uncertainty quantification is necessary. This 
is because the defect is defined in the highly localized regime and thus contains 
a high degree of intrinsic uncertainty. In our previous work (2l| . some linear 
wave equations with the singular source term were considered with the Legendre 
PC method. Our current work takes the first step to investigate the uncertainty 
quantification of nonlinear optics equations, with the source term confined in 
the highly localized regime. 

The paper is composed of the following sections. In Section 2, the Klein- 
Gordon equation with the singular source term is considered. In this section, 
we show that there exists a critical phenomenon and find the exact value of the 
critical parameter value. In Section 3, we briefly explain the gPC method and 
derive the gPC expansion of the Klein-Gordon equation using the Legendre poly- 
nomials. In Section 4, numerical methods used for the numerical approximation 
of the deterministic equations derived in Section 3 are explained including the 
consistent method for the approximation of the Dirac 5- function. In Section 5, 
numerical results of the gPC method for the Klein-Gordon equation are given. 
Comparison between the MC method and the gPC method is given. In Sec- 
tion 6, the critical phenomenon for the sine-Gordon equation is explained. The 
method of finding the critical initial velocity is explained using the PDF of the 
receiving signals and the gPC mean value. The Hermite gPC method is pre- 
sented. In Section 7, the Legendre gPC is used for the sine-Gordon equation 
and the deterministic equations are derived. Numerical results and remark on 
the error are presented. In Section 8, we briefly outlined the gPC method for 
the case that the randomess is in the amplitue of the singular potential function 
and explain that the same gPC method can be applied to this case. In Section 
9, concluding remark and future works are presented. 
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2. Critical phenomenon in linear equation with a singular potential 
term 



2.1. Inhomogeneous Klein- Gordon equation 

We first consider the one-dimensional, inhomogeneous, Klein-Gordon type 
equation in the general form given by 



where u,p 



utt - c^Uxx + d^u ^ p{x,t), a; e R, t £ M+, 
[ X M+ — > R, c, d e R. The initial conditions are 



(1) 



u{x, t) 
ut{x,t) 



u°{x), i = 0, 
0, t = 0, 



where vP ,u d !?■ 



u{x,t)->0, |x| -> -foo, t>0. 
The solution of Eq. ^ is given by 



u(x, t) 



1 



dT 



X-\~c{t — T) 



-c(t-r) 



[c^t-rf-ix-O' 



di, (2) 



where Jo is the Bessel function of the first kind. The above solution can be 
derived using the Green's function method 0. 

Now suppose that the source term involves the singular function such as 
p{^,t) = r]S{£,) where 77 S R is a constant and 5{^) is the Dirac ^-function. 
Then the solution becomes 



''-[c'it-rf-x^] 
c 



dr 



H{c(t 



(3) 



where H is the Heaviside function and c 7^ 0. As shown in the solution form, the 
parameter rj only scales the solution but does not yield the critical behavior of 
the solution. In order to obtain the critical solution, we consider the following 
several different cases. 



2.2. Steady-State Solution: c = l,d = l,p = r]S{x) 

We first consider the case of a possible steady-state solution with a singular 
source term. Consider the case that c — 1, d — l,p ~ riS{x) for the bounded 
domain x G Q = [—1,1] with the boundary conditions 

u{-l,t) = = ux{l,t). 

Then there exists a steady-state solution that satisfies 

-Uxx + u^r]S{x), w(-l,-)=0, Mx(l,-)=0. 
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By using these boundary conditions, the general solution can be easily derived 
as 



u{x) — 77 



1 



sinh(.T + 1) — H{x) sinh(a;) 



If different boundary conditions are chosen such as 
u(-l,-) = l, u,,(l)=0, 

then the solution becomes 



??e^(e^ + l) 
2 e4 + 1 



1 



1 



1 



2 e4 + 1 



These examples show that there exists the steady-state solution for any value 
of J] but we do not observe the critical phenomenon in the solution induced by 
the singular potential term p{x) = rid{x). 

2.3. Self interacting potential: c — l,d — 0,p{x, t) = V{x)u{x, t) 

For another possible steady-state solution, we consider the self-interacting 
potential term given by 



utt - Uxx = Vu, a; e ri, t > 0, 



(4) 



where : M ^ R is the potential function. We show that some V yield a critical 
phenomenon. 

Similar equation has been used in quantum physics. When V is given by the 
singular source term, the above equation represents some physical model such 
as one-dimensional diatomic molecule model with the singular potential well in 
the presence of static electric field 

The boundary conditions we choose are 



t) — 0, 

Ut +Ux = 0, 



X 
X 



-1, i>0, 
1, t>0. 



(5) 



with some non- vanishing initial condition 

u(x,0) = u"{x) ^ 0. 

The motivation of the given boundary condition of {dt + dx)u = at x = 1 is 
that if the potential term V is localized inside 51 and vanishes at the boundary, 
then the given boundary condition yields the outflow boundary condition. The 
potential term V will contain the singular term such that the potential effect 
vanishes at the domain boundary x S dVl = {±1}. 

To show the criticality of the solution, we define the energy E{t) of the 
system given by 

Eit) = \{^ ut\\l. + \\ Hi.). 
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The rate of change of energy E with respect to time t, E{t) is then simply 
computed by 

E{t) = j UtUttdx + [utUx\\ - j UtUxxdx. 

Using the boundary conditions we have 

E{t)^ J uutVdx - u'^t{l,t). 

Suppose that the potential function V is given by 

V ^f]S{x-a)u{x,t), r)>0, a e (-1, 1). (6) 

Then 



E{t,ri)=ri J uutS{x ~ a)dx — Uf{l,t) — riu{a,t)ut{a,t) — Uf{l,t). (7) 

Without loss of generality, we choose a = 0. Using the continuity of u at 
a: = a = 0, the jump condition is given by 

lim [u^(-e,t) -M^(+e,t)] = ?7u(0,t), Vi G IR+. (8) 
Lemma 1. The steady-state solution is determined by a unique rj'^ and rj'^ = 

Proof. Using the boundary conditions, we easily show that the steady-state 
solution u'^{x) is given by 

-'i-) = \ Sfl'} "^^7\f ^ (9) 

^ ' \ C{1 + a) a; e [a, 1] ' ^ ' 

where the constant C depends on the initial condition. Notice that the param- 
eter T] does not appear in the steady-state solution. Using the jump condition, 
we obtain 

e-i.O+ 

Using Eq. (|9]), the critical value of rj, rjc yielding the steady-state solution is 
obtained uniquely 

7?== hm [<(-e)-<(e)]K(a) = -^. □ 
e^0+ 1 + a 

Theorem 2. There exists a critical behavior around rj'^ obtained in Lemma 1, 
that is, 

E{t,r])>0 if7/>?7^ 
Eit,7^)^0 if77 = ^7^ 
E{t,'n)<0 if77<^7^ 

for t — > oo where rj G N^{r]'^) for an arbitrary small value of e > 0. 
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Proof. Without loss of generality we consider a = 0, i.e., the position of S- 
function is at the origin and suppose that we choose 

rj ~ jf + Ar/ = 1 + A77, 

where rj'^ = 1 (from Lemma 1) and < |Ar/| <C 1. Based on these assumptions, 
we look for the solution in terms of the power series in Ary, that is, 

u^u' + Aryui + iAr]fu2 + {Arifus H , 

where is the steady-state solution. Plug u into the PDE given by Eq. |3]with 
[6] and then each of (A?7)" terms yields 





{ui)tt 


- {Ui)xx 


= 5{x){u' - 


hUi), 


{^nf ■■ 


{U2)tt 


- (U2)xx 


= 5{x){ui - 








- (U'i)xx 




I-U2), 



We choose the initial condition such that the steady-state solution becomes 

1, for a; e [0,1] ' 

For the PDE of the first order term above, the boundary conditions and the 
continuity condition yields 



ui{x,t) 



a{l + x) + ct{l + x), forx<0 
a — cx + ct, for X > 



The jump condition, ~Ux{e^) +Ux{e~) = u*(0) + iii(0) yields c = 1. We choose 
the initial condition such that a — 1. Then by ignoring the higher-order terms 
(A??)("), n = 2,3,4, •• •, we obtain 



u{x, t) 



l + a; + A77(l + a; + f(l + a:)), for a; < 0, 
1 + Ar7(l - a; + t), for a; > 0. 



If A77 > the solution increases unboundedly and the energy increases with 
time. If A77 — 0, the solution becomes the steady-state solution and the en- 
ergy is conserved. If A77 < 0, the solution decays until t = — -I- 1^ For 

^ — ~ (sr; ^) ' higher-order term becomes significant and they should be 
considered so that the energy further decays. □ 

Remark 3. To get the complete solution, one needs to solve the higher-order 
terms. It may become difficult to solve those PDEs exactly but the proof is 
enough to show that there exists a critical behavior of the solution around rj'^ 
when \Arj\ <^ 1. For example, if Arj « 5x 10~^, the first- order term is dominant 
up to t 0.2 X 10^, which is enough for the long time behavior of the solution 
for our purpose. 
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Figure 1: Left: Analytical solutions of u{x,t) at different final times tf = 0, 500, 1000, 1500. 
Right: Contour plot of E for rj £ [0.9995, 1.0005] 



The energy of the system up to the first order term is given by 
And the rate of change of E with time is then given by 

E = viv-^)[v+tiv- 1)] + 1)' a+t). 

Figure [1] shows the solutions at different times tf = 0, 500, 1000, 1500 for 
different values of 77, rj — 0.9995, 1.0005 around rj'^ = 1 (left) and the contour 
of in 77 and t plane (right). 

2.3.1. Linearized sine- Gordon equation 

The same phenomenon can be obtained for the linearized sine-Gordon equa- 
tion with the point-like source term where sin(u) « u such that 

utt - u^x ^ ivH^) - '^h, a;e(-l,l), t>0, (10) 

with the same boundary conditions 

u{—l,t) — 0, and Ut + ^ at a; = 1. 

As in above, the energy, E(t) = ^ J^-^ (w^ + u^) dx, has the rate of change 
below 

E{t) = J -uut{x,t)dx + T]ut{0,t)u{0,t) - uf{l,t). (11) 

We have a similar phenomenon E{t) ^ J^-^ —uut{x,t)dx + r]u{0,t)ut{0,t) — 
it((l,i) = 0, that occurs around the following steady-state solution 

,,s(^^\C{l + e^){e-^-e-+^), for x < 0, , 

where C is a constant that depends on the initial condition. The critical value 
of 77 is Ty'^ = 2coth(2). 
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Table 1 : Correspondence between standard forms of continuous probability distributions and 
Askey scheme of continuous hypergeometric polynomials 24 , li3\ . 



Distribution 


Density Function 


Polynomial 


Weight function 


Support range 


Normal 


1 


Hermite Hn{x) 


e 2 


(— oo oo) 


Uniform 


1 

2 


Legendre Pn{x) 


1 


[-11] 


Beta 




Jacobi P„(a, /3)(a;) 




[-11] 


2i + a + f<_B(a+l,^+l) 


Exponential 




Laguerre Ln{x) 




[0 oo) 


Gamma 


r(Q + l) 


Generalized Laguerre l\^\x) 




[0 oo) 



3. Expansion using the gPC method 

3.1. gPC method and Askey scheme. 

We use the gPC method for the solution of the Klein-Gordon and sine- 
Gordon equations using the Wiener- Askey scheme [s^, in which Hermite, Leg- 
endre, Laguerre, Jacobi, generalized Laguerre orthogonal polynomials are used 
for modeling the effect of continuous random variables described by normal, 
uniform, exponential, beta and gamma probability distributions, respectively 
35|- These orthogonal polynomials are optimal for those PDFs since the 
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weight function in the inner product and its support range correspond to the 
probability density functions for those continuous distributions. 

Table (1) shows the set of polynomials which provide an optimal basis for 
different continuous PDFs [i^,!!^- It is derived from the family of hypergeomet- 
ric orthogonal polynomials known as the Askey scheme, for which the Hermite 
polynomials originally employed by Wiener [32j are a subset. 

Following the standard GPC expansion, we assume that u{x, t, ^) is suffi- 
ciently smooth in ^ and has a converging expansion of the form 



{u,x,0 = 2^Uk{x,t)Pk{0. 

k=0 

where the orthonormal polynomials Pk{0 correspond to the PDF of the random 
variable ^ and satisfy the following orthogonality relation: 



nPkPi] 



PkiOPiiOpim ^ Ski, 



where 6ki is the Kronecker delta and p(^) is the weight function. Note that the 
polynomials are normalized. 



3.2. Expansion using the Legendre gPC for the Klein- Gordon equation 

We use the gPC method to find the critical value of 77 numerically. Assume 
that ?7 > is a random variable and 77 G [a, h] has the uniform distribution. 
Let ^ be a random variable with the uniform PDF over £7 = [— 1, 1]. Let 77 be a 
linear map, 77 : [a, 6], 77 = -I- 
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The solution u{x,t,^) will be expanded using the Legendre polynomial in ^ 
with time t and the value of f will be sought that yields the invariant solution 
with time when t — > oo. Such ^ or 7y is then determined as the critical value for 
the Klein-Gordon equation. 

For this, we first expand the solution u in the Legendre polynomial in the 
Galerkin sense such as in |17| : 

oo 

u(a;,t,0 =5I"'(^'^)^'(0, (13) 

/=o 

where ui{x,t) are the expansion coefficients and the Legendre polynomials 
• We consider the truncated version of Eq. by considering the first iV + 1 
terms. Plugging Eq. ([T^ into Eq. Q and using V from Eq. ([5]) with a = 0, 
we get 

E [{"'(^' - {'^'(^' ^)}-] ^'(^) = -T^^ + ^ E ^'(^' ^)^'(^)- 

1=0 ^ i=0 

Multiplying both sides by yields 



JV 



1=0 



b — a ^ a + b 



N 



5{x)Y,ui{x,t)Li{i)Li^ 



1=0 



Integrating both sides over [—1, 1] and using the orthogonal property of the 
Legendre polynomials, we get, 



21 + 1 21 + 1 



5{x)ui{x,t) 



b — a 



1=0 



(14) 



For the second term in the right hand side of the above equation, one can show 
the following relation 



2(i+l) 

(2;+i)(2;+3) 

21 



(2;+i)(2;-i) 

otherwise 



for l^l +1 
for l = l' -1 



(for the derivation of the relation, see Appendix 2) 
Then, Eq. \1A\ can be rewritten as 



{ul{x,t)}^ 



{w/(a;,t)}^^ = K^) 
ui[x,t) + 



b ^ a 



21 ~ 1 



ui-i{x,t) 



b- a\ f l+l 
21 + 2, 



ui+i{x,t) 



(15) 
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Substitute Eq. into the left boundary condition of Eq. ^ to get, 



N 



1=0 

By using the orthogonal property of the Legendre polynomials, we have 

uii-l,t)^0, yi = 0,---,N. (16) 
Similarly the right boundary condition of Eq. ^ yields 

N 



1=0 

By using orthogonality we have 

{ui{ht)}^ + {ui{ht)},^0, V; = 0,---,7V. (17) 
We choose the initial condition as 

N 



(l + x) =^ii,(a:,0)L,(O. 



1=0 

Again using the orthogonal property, we get 

{l + x) J^^Li{Od^^^^ui{^,0)- 
This can be written as 

Mo(x,0) = (1 + a;), 

uiix,0) = 0, V/==0, •••,iV. (18) 

So from Eqs. (US]), and (HH]) we have {N + 1) equations for the gPC 

solution of the Klein-Gordon equation. 



4. Numerical methods for solving PDEs with singular source terms 

Consistent Chebyshev method 

The derived gPC equations for ui{x,t) in Section 3 will be solved numeri- 
cally. For the space derivative, the Chebyshev spectral collocation method is 
used based on the Gauss-Lobatto quadrature points Xj = — cos(j7r/m), j = 
0, 1, 2, • • • , m where (m -I- 1) is the total number of the quadrature points for 
Xj G [—1, 1]. As shown in Section 3 the derived gPC equations for ■ui{x,t) still 
contains the ^-function. Solving the PDE containing the singular source term 
properly is the crucial part for obtaining the accurate and stable numerical solu- 
tions for ui{x,t). To deal with this problem, the (5- function is approximated by 
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the consistent formulation proposed in [19| where the consistent formulation was 
proposed for the one-dimensional problems and it was shown that the method 
yields highly accurate results at the collocation points. It is not straightforward 
to extend such a formulation to higher dimensional problem, but as shown in 
the following sections and also shown in 2^ 21 , 22 1 , this formulation yields an 



accurate result for some PDEs and can be expanded to higher dimension by 
using the dimension-by-dimension approach. Moreover, the method is simple 
and easy to implement. The underline motivation is that the approximation 
of the (5-function is sought such that the steady-state solution can be obtained. 
This method is also referred as the well-balanced equation in the literature. The 
steady-state solution of the equation, in the previous section satisfies 

-ul,=ijSu^ix). (19) 

Let 5m be the approximation of S on the (m + 1) collocation points. Then 
the consistent formulation approximates 6m as the derivative of the Heaviside 
function H[x), using the same spectral derivative operator. If Hm denotes the 
Heaviside function defined at the collocation points Xj , then 5m is given by the 
following: 

5{x) = 5m := D£r,„, (20) 

where D g R('"+i) x is the Chebyshev derivative matrix and Hm. e M(™+i) ^ ^ 
given by Hm — (0, 0, • • • , 0, 1, 1, • • • , 1)^. To avoid any ambiguity, we choose the 
odd value of m such that for fc = (m — l)/2. 



Cfc < < ^fc+i, and {H„ 



0, for i < k 

1, for i>k + l 



That is, the location of the (5-function is off the collocation point. Then the 
consistent formulation of Eq. (|19l) becomes 



- D^C/ = r,AU, (21) 
where A is the (m -(- 1) x (m + 1) matrix whose diagonal elements are dm and 

U =iUixo)r--,U{Xm)V- 

4-2. Numerical errors of numerically obtained critical values of rj'^ 

As shown in Section 2, the critical value of r]'^ is determined by the following 

r?^= hm [<(-6)-<(e)]K(a). 

That is, the critical value is determined by the local values of u and near 
X = 0. Thus the convergence of the numerically computed value of rj'^ depends 
on the local convergence of u and Ux near x = 0. Figure [2] shows how the local 
convergence is obtained by our numerical approximation using the consistent 
method. In the figure, the convergence of the solution is obtained by solving Eq. 
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m+1 =32 




m+1 = 2048 



o 




Figure 2: Pointwise errors of \u — u"]. Left; Without the spectral filter. Right: With the 
spectral filter, m = 32, 64 ■ • • , 2048 from top to bottom. 



(PT|) iteratively. The left figure shows the pointwise errors with different values 
of m and the right figure the pointwise errors after the spectral filter is applied. 
As shown in these figures, the convergence can be enhanced away from the 
discontinuity when the filter operation is applied. But the overall convergence 
is slow due to the singularity at a: = 0, which is the Gibbs phenomenon for the 
spectral method This implies the slow convergence of the numerically 

obtained yy*^. 

To show such Gibbs phenomenon we consider the exact solution of Eq. (|19|) 
given by 

u-{x) = A[{x + l)-j^^H{Oddj. (22) 

The solution is continuous in (—1, 1) but at a; = it is not differentiable, u{x) € 
C°. Let the exact solution u'^{x) be expanded as the Gegenbauer polynomial of 
degree m: 

m 
fe=0 

where 

1 /"^ 

/ {l-x^Y-^C'^{x)u''{x)dx, 0<k<m, 
"fc J-i 

We use the particular case of the Gegenbauer polynomials for /i — 0, which 
corresponds to the Chebyshev polynomials. Figure [3] shows the Chebyshev 
Galerkin solution (left) and the Chebyshev collocation solution (right) with 
m + 1 = 128. As shown in these figures, there exist the Gibbs oscillations near 
X = 0. The Gibbs oscillation is seen mainly at the collocation point x — a;(m+i) /2 
for the collocation solution. These oscillations on the collocation points yield 
the error in determining the critical value rj'^. 

Table [2] shows the convergence of 77" with m, where 77" is the numerically 
determined critical value of 77*^. The critical value we will find from the following 
section using the given value of m will be those values in the table. For example, 
if we choose m — 63, the numerically computed critical value, 77" is rj" = 
1.004332 and our objective is to find that value as our critical value. That is. 
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exact solution 
+ numerical solution 



0.2 0.4 0.6 O.S 1 -0.5 0.5 1 

X X 



Figure 3: Left: The Chebyshev Galerkin solution. Right: The Chebyshev collocation solution. 



Table 2: Convergence of the critical value of rf^ with m. 



TO + 1 


77" 


32 


1.009595 


64 


1.004332 


128 


1.002036 


256 


1.000979 



given a finite value of to, tliere exists the corresponding 77" uniquely and when 
we find 77'^, we mean it by finding 77" instead of 77"^. 

4-.3. Time integration 

We use the second order centered difference in time for the time integration. 
For example let f/" e r(™+i)xi be the solution vector defined at t = n/S.t 
with the time-step 5t such as = {U{xo,nSt),U{xi,nSt), ■ ■ ■ ,U{xm,nSt))'^ 
for Utt — Uxx = P ■ Then the scheme is given by 

= 2[/" - [/"-I + (Jt^D • D[/" -1- P", 

where the element of the column vector P" e R(™+i)xi is (P")j = r]{Sm)jUj' 
and we apply D twice to U" for its second derivative in x. This is a two-step 
method in time. is given by the initial condition and we use the first-order 
approximation for U^. 

Then the numerical solution of ui G r(™+i)><i for each 1^0, 1, • • • , is 
obtained, by 

ul'+^ = [A{StyAiul'_^ + {BA{Stf + D\6tf + 21} ul' + A{5tf Ciu]^^^] - uf"^ 



where Ai = (^) P = Q = (^) and A = D • P,„. The 

boundary conditions are 

ur\-ht) = 0, 
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^ (i-AK(i,t) + Awr_i(i,t), 

where A = ^ and Sx = Xm — Xm-i- For the boundary condition at x = 1, we 
use the first order finite diff'erence both in time and space. 

5. Numerical results: determination of rjc and the mean solution 

5.1. Critical behavior 

For the numerical example, we consider (m+ 1) = 64 and the corresponding 
r] w 1.004332 from Table 2. 

In FigureH, we show E{t) versus time with r/ = {1.004332, 1.004330, 1.004334}. 
We find for rj = 1.004332 (blue), \E{t)\ < 0(10"*^) when t > 3000. This shows 
E{t) is almost constant with time. Another view of the steady-state solution is 
given in the right figure of Figure 21 where some wiggles are shown due to the 
inconsistency of the initial condition in the beginning but later the steady-state 
solution is obtained. Note that we use the initial condition u{x, 0, ^) = 1 + x. To 
verify the claim, we also calculate u{x, t) with r/ = 1.004332, 1.004300, 1.004364 
with different final time tf in Figure^ For r/ ^ 1.00432, tf = 1000, 3000, 5000, 
for rj = 1.004364 and 77 = 1.004330, tf = 1000, 2000, 3000. The figure clearly 
shows that the solution remains almost constant for r] — 1.004332. But for 
ri = 1.004300 the solution decreases monotonically and for 77 = 1.004364 the 
solution increases monotonically. 

5.2. Finding the mean value 

First we compute the mean value of the solutions with time. Let u(x,t) be 
the mean over ^. Then 

1 

uix,t)^-J ^uix,t,Od^- (24) 
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Figure 5: u(x,t) for r) = 1.004332 (left), r) = 1.004364 (middle), r) = 1.004330 (right). 



We consider three mean values obtained by the foUowing three methods. 

MC mean: To compute the above integral we use the uniform sampling of 
77 such that 

= Vo < m < ■ ■ ■ < Vi < ■ ■ ■ < Vm = 
and use the transformation ^ = ^_|_^'^^_ — . Then Eq. ([M)) takes the form 

u{x,t) = — u{x,t,^{T]))dri. 

We use the trapezoidal rule to calculate the integration. 

gPC mean: 

From Ea. ([M)) we have 



1 /-I ^ 

"(a;, = 2 / X! ^)^k{Cjd£. = ^ " t)8m = U(){x, t). 

•^"1 fc=0 

which is the first PC mode. 

Quadrature mean: We consider (iV+1) Gauss-Legendre (GL) quadrature 
points rji between [ri~ , ry"*"] and weights w^. To get mean value we have 

N 

We calculate the mean value of the solution u{x, t, rf) by the MC method at 
the final time tf = 100 with different uniform sample points 500, 1000, 2000, 
4000, 8000, 16000, 32000, 64000 for 77 e [0.95, 1.05]. Also we calculate the gPC 
mean at t/ = 100. The results are shown in Figured The MC mean solution 
converges to the gPC mean solution but the convergence is slow. As shown in 
the figure, more than 128000 MC samples are needed to match the gPC mean. 
This is because the energy growth/decay rate around the critical value of rj'^ is 
extremely asymmetric and the uniform sampling is obviously inefficient. 

The left figure of Figure [7] shows the convergence of the MC simulations. We 
calculated the mean solutions dXtf = 100 with the sample points Mj — 500 x 2^ 



1 2 
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where j — 0, • • • , 8. We plot ||i?j||2 versus Mj and Ej = PCavg — Uj. Here Uj 
is the mean solution of uniformly distributed Mj samples. The result is plotted 
in logarithmic scale and the decay is linear. 

The right figure of Figure [7] shows the mean solutions with the Legendre 
gPC, MC simulation with 128000 equally distributed sample points and the 
quadrature mean with 50 Gauss-Legendre quadrature points in [0.9,1]. The 
solid red line represents the result by the MC simulation which deviates from 
the gPC mean or GL quadrature mean. The GL quadrature mean with the small 
number of points coincides with PC mean. Using the quadrature sampling for 
the mean solution was similarly used in 




5. 3. Finding critical value of rj 

We first define a solution sequence U'^ {rj) at x = 1. The index k denotes the 
sequence corresponding to time tk such as 

k:tk^ kAt, At > 0. 

Define ry', I — 1, 2, 3, • • •, 

v'^{v\u\v)-u^-\v) = o}- 
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Figure 7: Left: Errors of the MC mean with different number of sample points M, where 
M = 500 X 2^ for j = 0, 1, 2, ■ ■ ■ , 8. Right: Comparison of u{x, 100) derived from the MC 
simulation with 128000 sample points, the Gauss-Legendre quadrature with 50 quadrature 
points (red circles) and the Legendre chaos method. It is observed that the quadrature mean 
and the PC mean are almost coincident. 

For large t, U^{rj) is a monotone function with respect to 77 around the critical 
value of 77. As there exists the critical value of 77, rjc , there exists 77 such that 
U^{ri) — U^^^{ri) = 0. Now we define another sequence At; such as 

A?7' =?7'-77'-i. 

As [/' is monotone with respect to rj for large t, we have 

hm Ar/ = 0. 

and 

lim r/' — rjc- 

1^00 

Figure [8] shows the sequence of C/' and ry' by the MC simulations with 10^ 
sample points at different time intervals. The final time varies from 50 to 600 
and we calculate the solution at tf = 50, fOO, 200, • • •, 600. For any fixed tf, 
we plot the solution for different rj. If 3 77c then the graph u(l,tf,ri) converges 
to zero as t — > cx) when rj < rjc and similarly ty, ry) diverges to 00 as i — cx) 
when > rjc- From Figure |S1 we find that the solutions intersect at the common 
point {r]^,Us{l,tf,ril^)) and we claim that such 7;" is the critical value of rj and 
Us{l,tf,r]^) is the steady-state solution. 

The bottom figure shows the magnified region near 77" . If we see the figure 
carefully, it is then observed that the solutions for = 400, tf = 500 and tf = 
600 intersect almost at the same point, whose coordinate is approximately equal 
to (io" 0"i87728 ;^qo.30io295) ^ (1.004331952841278,1.999997717384314). Then 
our estimates are 77;? = 1.004331952841278 and ti^(l, t/, 77^) = 1.999997717384314, 
In Section 4 we found for (771 + 1) = 64, the critical value of 77 is about 1.004332, 
which is close to the value of 77". Here note that the value of 77^ in Table 2 was ob- 
tained by many of the individual simulations of the equation with different value 
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Figure 8: Top: The locus of u{\,tf^r\) with different final time by the MC simulations, r] 
varies from 0.95 to 1.05. The curves seem to intersect at the common point. The lowest curve 
(blue solid line) in the right side of intersection point is for tj- = 50 and then the top curve 
(black solid) for tf = 600. Bottom: The close view around the intersection point. For the 
MC simulations, the intersection points are distributed randomly in a very small region but 
for higher time the curves tend to converge to a single point. 



of ?7, which is highly expensive and rf was obtained with only 10 6 accuracy. 
The error for the steady-state solution at a; = 1 is 2 - 100-3010295 ^ 2.2 x 10"^, 
which agrees our error estimation in Section 5.1 and FigureU) The error of 10~^ 
conies from the Gibbs osciUations at a; = 1 as explained. 

Figures [S] shows how to find the critical value rjc by the gPC method. 
In the gPC method we reconstruct the solutions by considering the first 80 
modes, at different time t/, where tf runs from 50 to 350 with the interval 
length 50. At every stage we plot u(l,t/,^), where ^ G [—1,1] and con- 
sider 10^ equally distributed sample points generated by the 80 gPC modes. 
We found that the solutions intersect at the common point (^c, i/, fc)) 
and we estimate that it is the critical value of ^ for S [0.95, 1.05]. We 
find the approximate value of is = 10-i"62285 _ 8.6639313 x 10"^ and 
Ws(l,*/,?c) = 10°-30i032 ^ 2.000009230329. For this value of the corre- 
sponding rjc is rjc — 1.004331965650747. The value of r/^ divides the interval 
[0.95, 1.05] into two subintervals [0.95, ?7c] and [?7c,1.05]. When i] < rjc, the 
graph converges to zero as tf — )■ oo and similarly the graph diverges to oo if 
r] > rjc sls tf ^ oo. In Figure [9l the right figure is the magnified version of 
the left figure near ^c- The graphs intersect within a very small region between 
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-0.05 0.05 0.1 0.15 0.2 0.25 jQ-i.uf.^J ^^-im^j 

Figure 9: Left: The locus of u{l,tf,ri) with different final time by the Legendre chaos, r] 
varies from 0.95 to 1.05. The curves intersect at the common point. The lowest curve (cyan 
solid line) in the right side of intersection point is for t y = 50 and the top curve (black solid) 
for tf = 500. Right: A close view around the intersection point. For the Legendre chaos 
method, the intersection points are distributed randomly in a very small region. Note that 
the horizontal axis is for ^ not -q. 



jqO.301028 ^ 1.999990809648848 and 100-30i032 ^ 2.000009230329776 in the ver- 
tical axis and between 10-106229 ^ 8.663831554873760 x lO'^ and IO-106228 ^ 
8.664031049264383 x 10"^ in the horizontal axis. The last two lines aXtf = 300 
and = 350 (blue and black) intersect at Us = 10°-3°i°32 ^ 2.000009230329776 
and the error in this case « 7 x 10~^. The error in calculating the critical value 
of r] must be less than the difference of io~i 06229 -[^q-i.06228 -^^j^ich is almost 
equal to 1.99 x 10~^. Thus we get the result with in the accuracy of ^ 10^^. 
The difference between 77^ computed by the MC method and the gPC method 
is 

At^^ \r]''{gPC) - r]''{MC)\ - 1.2809 x 10"^ 

Here note that for the gPC method we only use = 80 while for the MC 
method we use M = 10^, which shows that the gPC method can determine 77'^ 
efficiently. 

6. Critical phenomenon of the sine-Gordon equation 

6.1. Basic idea and problem formulation 
Consider the sine-Gordon equation 

Uxx ^ + sin(M) =0, X e M, 

where u : K x M+ — > M. The soliton (or kink) solution of the sine-Gordon 
equation is given by 



/ N . -if f 

u(x,t) =4 tan < exp , > 
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This represents the continuous profile with u as a: — ;> — oo, and u — !> 27r 
as a; — +00. This solution propagates with the constant velocity V. Another 
solution can be obtained in the form 



u{x, t) ~ A tan 



exp 



Vt 



-11 f x-Vt 
= 4 cot < exp 



which is known as the antisoliton(or antikink). 

We consider the sine-Gordon equation with an impurity 



Utt — + sin(u) = e5(x) sin(M), x £ 



t £ 



(25) 



where e € is the amplitude or strength of the singular potential function. 

The full sine- Gordon equation with an impurity Eq. (j25p is, however, non- 
integrable and no exact solution is known. The problem we solve numerically 
in this work is defined by Eq. ([25]) in x £ il, where il — [-~L,L], L > with 
the initial kink velocity V. For simplicity we assume that the singular potential 



exists at cc = 0. For the numerical computation, we set L = 8, e 



and the 



initial kink front xq = —6. The initial and boundary conditions are given below. 
Initial conditions (kink soliton) 



u(x, 0) 
utix,0) 
Boundary conditions: 



4 tan 



-1 



W 



f X- Xq 

,i + -p(?»). 



-V2 



ut{L,t) + Ux{L,x) = 0, 

u{~LA) = 4tan~ 



exp 



-L- Xq- vt 



The above boundary condition is used for the outflow boundary condition. 
Note that this boundary condition is not an exact outflow boundary condition 
because of the nonlinear potential term, sin(u). Better boundary conditions and 
exact outflow boundary conditions can be found in [i^l • Since the main goal of 
our study is to find the critical velocity, the above boundary condition does not 
affect the results. 

Figure [To] shows the solution behaviors around the critical value of Vc, 0.1 < 
Vc < 0.2. The left figure shows the sub-critical solution behavior when V = 
0.1 < Vc and the right figure when V = 0.2 > Vc- The horizontal axis denotes 
the time and the vertical axis the location of the kink front (the red color 
shows the kink front and the 2tt plateau). These solutions are obtained by the 
direct numerical approximation using the Chebyshev collocation method with 
the consistent approximation of the 5- function explained in the previous section 
with m + 1 — 128. As shown in the figure, the soliton solution is captured in 
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Figure 10: Spectral solution of the sine-Gordon equation with m = 127. Left: Sub-critical 
behavior with V = 0.1. Right: Super-critical behavior V = 0.2 [20| . 

the potential well ii V < Vc while it eventually passes through the potential 
well. That is, at a; = L, u{x,t) becomes either 27r (particle-capture or particle- 
reflection) or (particle-pass) as i — >■ oo. We use this solution behavior to 
find the critical velocity Vc- Here note that the early numerical study of the 
discretized sine-Gordon equation is found in the literature (28[. 

6.2. Finding the critical velocity 

In the previous section, we considered the Klein-Gordon type linear equation 
for which we use the MC sampling using the gPC expansion for constructing 
the solutions. However, for the nonlinear sine-Gordon equation, such a MC 
sampling is not usable because of the two extreme PDF of the solution for large 
t. The construction of the gPC for large t yields either trapped or transmitted 
solution. 



(1 - Vc)S{t) 



VcS{t - 2tt) 



2tv 

z = u°° 

Diagram 1: The PDF of the solution for t —J- cxd, at x = L for uniform V. 



The critical value can be computed by the weighted ratio of the two extreme 
PDF as shown in Diagram 1 if we know the PDF of V. The ratio is indeed the 
amount of how much the mean value deviates from the value of 2t: for the trapped 
solution. Thus the critical value of V can be still found from the gPC expansion 
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using the mean mode. If F € [Va,Vh] and if Vc is the critical velocity, then for 
V £ [Va,Vc], the solution is 2tt (because the solution is trapped), and for V S 
(VcH], the solution vanishes (because the solution is passed). If f {V;Va,Vb) 
is the PDF for V, then 

u{x = L,tf) = / f{V-Va,Vb)u{x,tf,V)dV 

= 2^ / f(V-Va,Vb)dV. (26) 

For example, if we use the Legendre chaos, the distribution function is uniform 
and f {V;Va,Vb) = y^^y ■ For this by using Eq. (^5]). the Legendre mean uieg 
is given by 

2tt 

Uleg = — {Vc - Va) , 

Vb — Va 

and we have 

Vc = Va + {Vb-Va)^. (27) 

6.3. Hermite chaos and the distribution function of input velocity V 

Suppose that the input velocity, V, is the random variable which is normally 
distributed with mean /i and standard deviation a. Here V, conditional on 
a < V < (3, has a truncated normal distribution. We consider ^ = ^ {a + ^) 
and a ^ 10^^ where a — Va, P — Vb. The reason for choosing small a is to keep 
the error - 10"^. The PDF of V, f, is given by 

' > for a<V <P, 



f{V■,^l,a,a,P) = l WF^(^' "-"-A^' (28) 

[ 0, otherwise. 

The density function involves '0(-) which is the PDF of the standard normal 
distribution and $(-)j its cumulative distribution function 4]. Due to the small 
value of a, 

' e-'-^H, ( ^) Hr ( ^) dV « aV2^U S,, , 



for which see Appendix 1. 

6.4. Sine-Gordon equation and Hermite chaos 

We expand u{x,t, V) in terms of the Hermite polynomials as follows 



,v-i^. 

ui{x, i)ni\ 

1=0 



iix,t,V) = ^M;(a;,i)^^;(- 
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Figure 11: Left: The mean solution and the standard deviation at t final = 600. Right: The 
first few PC modes for the Hermite chaos at tj^^^i = 600 and I = 1, ■ ■ ■ ,7 . 



For the numerical purpose, we consider the first + 1 modes. Then the sine- 
Gordon equation can be written as 



[ui{x,t)\tt 
where ipiixTt) 



t! V 2,71 

N 



and 



r sin [Y,ui{x,t)Hi[^)\Hi{^)^{-i)d^, 



For the initial condition, we have 

f X- Xq 



4 tan 



exp 



^ui{x,^)Hi 



1=0 



By using the orthogonal property of the Hermite polynomials we get 



ui{x, 0) 
Similarly 



V2^U J- 



-4 



tan 



exp 



X — Xq 



Hi{j)uj{'-f)dj. 



27r/! 



{ji + CT7) 



exp( , J 



1 - (/u + (77) 1 + exp 



2(x-xo) 



■Hi{'y)uj{-/)d'y. 
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The boundary conditions are given by 
{ui{L,t)}t + {ui{L,t)}^=0, 



ui{-L,t) 



for 1 = 0,1, ■■■,N. 



27r/! 



tan 



exp 



-L - xo - (/u + 0-7) t 



1 - (m + o'T) 



Hi{"f)uj{j)d'y, 



6.5. Mean and Standard Deviation of the solution 
By definition, the mean solution is given by 



u{x,t,V\VG[a,/3]) 



f{V;a,(3,n,a)u{x,t,V)dV 
1 



2Tra 



$ (^^^ - i> (^) 



e 



y^jii{x,t)Hi 



dV 



N 



uo{x,t) + 



uo{x,t) + 



$ (^) - $ i^)] 



1=1 

N 



^ui{x,t) e ^2 if; (7)^7 



$ (^) - $ i^)] 



By using the orthogonaUty of Hi (7) 



To calculate standard deviation we first calculate u'^{x,t). 
u^{x,t) = / fiV;a,0,ii,a)u\x,t,V)dV 

J a 



(T 



fiV■,a,P,^I,a)ulix,t)H^{—^]dV+ / f {V;a, ii,a) ul{x,t)Hl 



V ■ 



+ ■■■ + '^Y.Y.I f{V\a,p,ii,a)um{x,t)un{x,t)HmO^^]Hr.\^^\dV 



AT 



^ m\ul,{x,t), 
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and 



[V.ulix, t) + 2\ul{x, t) + ---+ Nlu%{x, t)] 



The standard deviation is then given by 
a = \J u^{x, t) — {u{x, t)}'^ s 



AT 

E 

m— 1 



mlul,{x,t) 



6. 6. Estimation of critical velocity 

To estimate the critical velocity we use the Eq. (PE)) . Here the density 
function / is given by Eq. (|28p . By using these two equations we get, 



Uhc 

2tt 



1 



a — ji 



(29) 



If Uhcr = 27r, $ 

trapped and if 



0, $ 



= $ ^^^^ ) which implies Vc — (3, all solutions are 



$ (^^^) which implies Vc = a, all solutions 
are passed. For numerical simulations, we consider /x = 0.12, ct = 0.01, a = 
0.11, (3 = 0.13. We obtain Uhor « 3.24341621 at the final time tf = 600. By 
using Eq. ((29)) . we get. 



Vc-0.12 
0.01 



3.24341621 
2^ 



0.13-0.12 
(101 



0.11 - 0.12 
OOl 



0.11-0.12 
OOl 



Based on these results and by using the statistical table [3] we estimate Vc ~ 
0.1203. This estimate is close to what we get from the MC simulations and 
the Legendre chaos in the next section. The obtained value is, however, not 
accurate; it matches with the value obtained by many MC simulations within 
two digits only. Figure [TT] shows the initial condition, the mean solution, and 
the standard deviation in the left figure and the first 7 PC modes {I = 1, ■ ■ ■ ,7) 
in the right figure. 



7. Sine-Gordon equation with the Legendre chaos 

Suppose that V G [Va, Vb] = [a, b] has the uniform distribution and 

= + ee[-i,l]. 

We consider u{x,t,£^) — '^f^oUi{x,t)Li{£^). By using the orthogonal property 
of the Legendre polynomials and the initial and boundary conditions we have 
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ui{x, 0) — 2{2l + 1) J^^ arctan |exp 
[u,(x,0)], = -2(2? + l)/^,^Z|L= 

ui{—L, t) = 2{2l + 1) J^^ arctan l^exp 
{uiiL,t)}^ + {uiiL,t)}^ = 0, 



f oxp 



i+exp( 

-L-a;o-tV(g) 



for ? = 0, 1, • • • , iV. Then we get 

N 



N 



E [{^'(^' - ^)}xJ ^Ke) = US{x) - 1} sin ^ i)ii(e) 



i=0 



for I — 0,1, ■ ■ ■ , N. Using the orthogonal property of the Legendre polynomials, 
we get 



{Mx,t)}tt-'{ui{x,t)}^^ = 
for Z = 0,1, •••,7V. Or 

{Mx,t)}tt - {Mx,t)} 



{e5{x) - 1} I'^sin [Ylui{x,t)Li{0^ Li{(,)d{0. 



21 + 1 



{e5(x)~l}Ux,t), 



for ? - 0, 1, • • • , TV and t) = j\ sin (j^^^^ ui{x, t)Li{0) Li{£,)d{0. 
The discretized equations are given by 



-li+i 



D' (^r) + (^)(e(D^i/)-/)0r 



[uz(x,0)]j 

ui{-L,t) 
ur\L,t) 



-2{2l + l) 



exp 



L-xo~ tViO 



2(2Z + 1) J arctan < exp 



r(i,0- (^) [u?{L,t)-u-+\L~-5x,t)] 



for I = 0,1,---,N. Figure shows the PC simulation with V^^ = 0.1215 and 
Vfc = 0.121757. The left figure shows the mean solution and the standard 
deviation and the right figure the first f 5 PC modes. The first mean mode is 
plotted with the red circle. From these results we find the mean solution uieg, 
at X = i. uieg = 2.011976945534794 at i/ = 600. By using Eq. we 
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Figure 12: Left: Mean and Standard deviation of the gPC solution at the final time, tf = 600. 
Right: First 15 gPC modes. Red circles represent the 1st PC mode (mean solution) 
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Figure 13: Solutions at z = 8 with time for V = 0.121582295531 (sub-critical) and V 
0.121582295532 (super-critical). 



obtain Vc ~ 0.1215822955316, which is close to the critical value by the MC 
simulation, which was 0.121582 < Vc < 0.121583 To verify the obtained 

result, two direct spectral simulations are conducted with Vi — 0.121582295531 
and V2 = 0.121582295532. Figure [T3] shows that the sub-critical solution is 
obtained for Vi ~ Vc — 5.9999 x lO^^'^ and the super-critical solution for V2 = 
Vc + 4.0001 X 10^^"^. This verifies that the obtained Vc is as accurate as at least 
up to 13 digits. 

Remark 4. For the Hermite chaos, the integral for the mean solution is not 
exact and we need to rely on the statistical table to get the probability values. 
For the Legendre chaos, however, we could find a highly accurate result. The ar- 
tificial oscillations fluctuating aroundO (particle-pass) and2n (particle- capture) 
shown in the individual solution at x — L due to the non-exact boundary con- 
dition and the Gibbs phenomenon are highly reduced in the mean solution and 
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the analysis shown in Diagram 1 yields an accurate mean solution resulting in 
finding an accurate critical velocity Vc- 



8. sine-Gordon equation with the uncertainty in e 

Suppose that we have the uncertainty in e, and e e [a, h] has the uniform 
distribution. Then 



b — a a + b 



We consider 



N 



u{x,t,^) = ^ui{x,t)Li{^). 



1=0 



uo{x, 0) = 4 arctan |exp 
Um{x,0) = 0, 



X — Xq 



[uo(a;,0)] 



cxp . 

iV ) V Vi-V' 



[u„(a;, 0)]j = 0, 
UQ{—L,t) =4 arctan |exp 
Um{-L,t) = 0, 



-L-Xn-tV 



]} 



{^;(i,i)}t + {«/(i,i)L = o, 

In the similar way as in the previous section we have for m = I,-- - ,N. 



N 



X [{Um{x, - {Um{x, t)}^J i„(0 = | ( + ^ 

TO=0 ^ ' ■' \m=0 

Using orthogonal property of the Legendre polynomials again, we get 



6{x) - 1 [ sin I ^ Umix,t)Lm{0 



{um{x,t)}^^-{um{x,t)}^^ 
for m = 0,1, - ■ ■ ,N. where 

<l>m{x 



2m+l\ (a + b^, . 1 , 6-a, 



■>Pm{x,t) 



,t) = j ^sm(j2ui{x,t)Li{0^L^{^)d{0, 



Mx,t)Li{0]L^iOd{0- 



The discretized equations are given by, 
form = 0, l,---,Ar. 
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uo{—L,t) = 4arctan|exp 
u,n{-L,t) = 0, 

u^^+^m) = u^{L,t) - (H) - u^+\L-5x,t)] , 
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Figure 14: MC simulations for different e to get the idea of the position of Cc. 

8.1. Finding critical amplitude 

To get the critical value of e, Ec, we follow the same procedure as in section 

6.2. Since e is uniformly distributed in [e(j,eti], we use Eq. (j27p . but in little 
different way. Let ea < ^ e& smd for e e [e^, ec], the solutions are passed and 
for e S (ec, eb], the solution get trapped. So Eq. (|27p takes the form 

Figure ([T^ shows some solutions at a; = with time for different values 
of e, e = 0.4925, 0.4975. For the captured solution, the solution at x = is 
highly oscillatory but for the transmitted one it is close to zero. From the figure 
we know that 0.4925 < Vc < 0.4975. Suppose we know for = 0.495 solution 
passes through the barrier and for £{, = 0.4975 solution is trapped by the barrier. 
By the gPC simulations fFigure[T5|) for = 600 we have uieg = 2.096586 . Using 
these in Eq. ^ we get = 0.4958342050637932. 



'L-XQ-tV 
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Figure 15: Mean and Standard deviation of the PC solution at Tf = 600 where e is stochastic 
variable. 

9. Conclusion 

In this paper, we proposed an efficient method using the gPC method to 
find the critical parameter values and the related statistical quantities for the 
Klein-Gordon and sine-Gordon equations with a point-like singular potential 
function. By assuming the unknown critical parameter values as a random 
variable, and expanding the solution in the orthogonal polynomials associated 
with the random variable, we computed the critical parameter with much less 
computation than the MC approach. For the Klein-Gordon equation, the critical 
parameter is associated with the amplitude of the singular potential function. 
Using the gPC reconstruction of the solution, the invariant point is found to be 
the critical parameter value. For the sine-Gordon equation, the initial velocity 
of the soliton solution is used as the random variable. For this case, both the 
Hermite and Legendre gPC methods are used. Due to the nonlinearity and 
the (5-function like PDF of the receiving signals, the reconstruction could not 
be used; the gPC mean was used instead. The Legrendre gPC mean converges 
quickly, and it efficiently and accurately determines the critical parameter value. 

As our numerical results show, the gPC method offers several advantages in 
finding the critical value as compared to the MC method. First of all, one can 
avoid the large number of individual simulations for different parametric values. 
Our proposed method suggests that, with a small number of MC simulations, 
it is possible to make the first guess about the critical value. Once the first 
guess is made, one can apply the gPC method to narrow down or pinpoint the 
critical value accurately and efficiently. Moreover, regarding the gPC method 
for the sine-Gordon equation, the stochastic variable appears only in the initial 
and boundary conditions, and does not appear in the main PDE. To perform 
the integration related to V, one can use the Legendre-Gauss quadrature for 
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the Legendre chaos and Gauss-Hermite quadrature for the Hermite chaos and 
obtain a fast convergence. In this paper, only 30 quadrature points are used. 
Adopting the quadrature rule, one can easily speed up the whole process. 

During our study, we also found a limitation of the gPC method regarding 
the determination of the initial critical velocity or the critical amplitude of the 
singular potential term for the sine-Gordon equation. That is, we could not 
use the reconstruction constructed by the gPC modes to find the PDF of the 
receiving solutions due to the non-linearity and non-separability related to the 
Gibbs phenomenon of the spectral method. 

Despite this limitation, we find that the gPC method works well for finding 
the critical values for the singularly perturbed Klein-Gordon and sine-Gordon 
equations. Since no significant gPC analysis exists for the singularly perturbed 
PDEs by the Dirac (5-function, our study could be a good resource for this type 
of research. In our future work, we will apply a similar gPC method for differ- 
ent PDEs, such as the nonlinear Schrodinger equations with singular potential 
terms. We will also consider more general types of uncertainties associated 
with the singular potential term, for which more than one random variables are 
involved and the gPC expansion should be in the multi-dimensional random 
space. To reduce the complexity due to the dimensionality, we will att emp t to 
use the gPC collocation method with the fast sparse grid methods 2^ 34, 36j . 
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Emmanuel Lorin for their useful discussions and references. 

Appendices 

Appendix 1: Hermite polynomials 

Hn{x) are the Hermite polynomial for ti = 0, 1, 2, 3 • • • and they are orthog- 
onal with respect to the weight function w{x) = e^~; 



H„iix)H„{x)e ^ dx — v27m!(5„j„. 

First few Hermite polynomials are Hq{x) = l,Hi{x) ~ x,H2{x) — x'^ — 
l,H3{x) = - 3x. 

The following few recursion relations are used in the paper: 

Hn+i{x) = xiJ„(a;) - 
H„{x) = nHn-i{x). 

71 

Rnix + y) - ^C,"a;'=i/„_fc(2/). 

fe=0 

Appendix 2: Legendre Polynomials 

Some important properties of the Legendre polynomials are used in the pa- 
per. Legendre polynomials are defined by the solutions Pn{x) of the following 
Sturm-Liouville equation: 

^[{l-x')^]+n{n + l)p,, = 0, xe[-l,l], (.1) 
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where Pn{x) is the Lcgcndrc polynomial of degree n. First few Legendre poly- 
nomials are po{x) = 1, pi{x) = x, P2{x) = ^(Sa;^ — 1), P3{x) = \{^x^ — 3x), 
Pi{x) = i(35a;'^ — SOa;^ + 3). Also = 1, 1) = 1 if n is even and = —1 
if n is odd, also p„(0) = if n is odd with the orthogonality 



Pm{x)pn{x)dx = ^ ^ ^mn- 



The recurrence relation is given by 

(n + l)pn+i{x) = (2n + l)xpn{x) - npn-i{x), (.2) 

P'n+lix) - P'n-lix) = (2n + l)Pn{x). (.3) 

From Eq. (.2) we have, 

(2n + l)xPn{x) = nP„-i{x) + {n + l)P„+i{x). 
Multiplying both sides by P^ and integrating from —1 to 1, we get, 

j xPn{x)Pm{x)dx = ^ j Pn-l{x)Pm{x)dx+^^^ j Pn+l{x)Pm{x)dx. 

By using the orthogonality of the Legendre polynomials we get. 



1 

xPn{x)Pm{x)dx 



2(m+l) 



p if n = m + 1, 



(2m+l)(2m+3) ' '<- — "■' i 

(2^+iK2m-i) ^ if n = m-l, (.4) 
0, otherwise. 



By using Eq. (.3) we get the foUowings. 
For n is odd number: 

p^+i = (2n + l)pn{x) + (2n - 3)p„_2(a;) + • • • + 7pi{x) + 3pi(a;). (.5) 

For n is even number: 

p^+i = (2n + l)pn{x) + (2n - 'i)pn-2{x) + ■■■ + 9pa{x) + bp2{x) + po{x). (.6) 

The derivative of Legendre polynomials as the finite linear sum of Legendre 
polynomials are used in the paper and some examples are as follows: For n = 6, 

P7(a;) = 13pe{x) + 9p4{x) + 5p2{x) +po{x). 

And for n = 7, 

Pg{x) = 15p7{x) + llp5{x) + lps{x) + 3}Ji(a;). 

By orthogonal property P'j{x)p^{x)dx = 0. So in general if m is even and n 
is odd or n is even and m is odd, 



1 

Pm{x)Pn{x)dx = 0. (.7) 

1 
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If both m and n are even numbers and m > n then by using the orthogonal 
property of Legendre polynomials and from Eq. (f75|) we have, 



1 



2(n-l) + l ' ' 2{n-3) + l 



7^— ^— + 3^- ^ 



2-3 + 1 2-1 + 1 

= 2[3 + 7+--- + (2n-l)] =n(n+l). 

If both m and n are odd numbers and m > n then the orthogonal property and 
Eq. CS]) yield, 

/ p^{x)p^{x)dx = 2[l + 5 + 9 + --- + {2n-l)]=n{n + l). (.8) 
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